/*
 * Copyright (c) 2025 CAX Conversion Project
 */
package com.comac.ins.isfd.util;

import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;

import java.util.Arrays;

/**
 * 铺层等效材料计算
 *
 * @author liji
 * @date 2024-8-8
 */
public class LayupCalculateUtils {
    public static void calculateEquivalentMaterial(String plyCode,
                                                   double singleThickness,
                                                   int symmetrical,
                                                   String layupSequence,
                                                   double[] compositeMaterialProperties,
                                                   double[] equivalentProperties
    ){
        // 计算特殊结构：N-蜂窝层，M-多层角片
        if (plyCode.contains("N")) {
            return;
        }
        else if (plyCode.contains("M")) {
            // 转化为铺层序列数组
            String[] layupStrings = layupSequence.split(",");
            int j = 0;
            double[] E1m = new double[layupStrings.length];
            double[] E2m = new double[layupStrings.length];
            double[] Nu12m = new double[layupStrings.length];
            double[] G12m = new double[layupStrings.length];
            double[] tm = new double[layupStrings.length];
            double[] Rhom = new double[layupStrings.length];
            double[] alp1m = new double[layupStrings.length];
            double[] alp2m = new double[layupStrings.length];

            double[] AE1 = new double[layupStrings.length];
            double[] AE2 = new double[layupStrings.length];
            double[] AG12 = new double[layupStrings.length];
            double[] ARho = new double[layupStrings.length];
            double[] AAlp1 = new double[layupStrings.length];
            double[] AAlp2 = new double[layupStrings.length];

            for (int i = 0; i < layupStrings.length; i++) {
                String layup = layupStrings[i];
                String[] paramList = layup.split(",");
                if (paramList.length == 2) {
                    int skinNum = Integer.parseInt(paramList[1]);
                    // todo equivalentList
                    double[] equivalentList = new double[8];
                    E1m[j] = equivalentList[0];
                    E2m[j] = equivalentList[1];
                    Nu12m[j] = equivalentList[2];
                    G12m[j] = equivalentList[3];
                    tm[j] = equivalentList[4];
                    Rhom[j] = equivalentList[5];
                    alp1m[j] = equivalentList[6];
                    alp2m[j] = equivalentList[7];
                }
                else if (paramList.length == 3) {
                    int metalNum = Integer.parseInt(paramList[1]);
                    // todo metalMaterialList
                    double[] metalMaterialProperties = new double[8];
                    E1m[j] = metalMaterialProperties[0];
                    E2m[j] = metalMaterialProperties[1];
                    Nu12m[j] = metalMaterialProperties[2];
                    G12m[j] = metalMaterialProperties[3];
                    tm[j] = metalMaterialProperties[4];
                    Rhom[j] = metalMaterialProperties[5];
                    alp1m[j] = metalMaterialProperties[6];
                    alp2m[j] = metalMaterialProperties[7];
                }
                AE1[j] = E1m[j] * tm[j];
                AE2[j] = E2m[j] * tm[j];
                AG12[j] = G12m[j] * alp1m[j];
                ARho[j] = Rhom[j] * alp2m[j];
                AAlp1[j] = alp1m[j] * tm[j] * AE1[j];
                AAlp2[j] = alp2m[j] * tm[j] * AE2[j];
            }

            double thickness = Arrays.stream(tm).sum();
            double ExxEquivalent = Arrays.stream(AE1).sum() / thickness;
            double EyyEquivalent = Arrays.stream(AE2).sum() / thickness;
            double GxyEquivalent = Arrays.stream(AG12).sum() / thickness;
            double Rho = Arrays.stream(ARho).sum() / thickness;
            double alp1Equivalent = Arrays.stream(AAlp1).sum() / (thickness * ExxEquivalent);
            double alp2Equivalent = Arrays.stream(AAlp2).sum() / (thickness * EyyEquivalent);

            // 保存信息
            equivalentProperties[0] = ExxEquivalent;
            equivalentProperties[1] = EyyEquivalent;
            equivalentProperties[2] = GxyEquivalent;
            equivalentProperties[3] = thickness;
            equivalentProperties[4] = Rho;
            equivalentProperties[5] = alp1Equivalent;
            equivalentProperties[6] = alp2Equivalent;

        }
        // 其他一般情况
        else {
            // 转化为铺层序列数组
            String[] layupStrings = layupSequence.split(",");
            // 如果对称，拓展数组
            String[] layupStringArray = expandSymmetrically(layupStrings, symmetrical);

            int[] layupIntArray = new int[layupStringArray.length + 1];
            double[] layerThicknessArray = new double[layupStringArray.length + 2];
            for (int i = 0; i < layupStringArray.length; i++) {
                layupIntArray[i+1] = Integer.parseInt(layupStringArray[i]);
                layerThicknessArray[i+1] = singleThickness;
            }



            double E1 = compositeMaterialProperties[0]; // 单向带沿纤维方向弹性模量
            double E2 = compositeMaterialProperties[1]; // 单向带垂直纤维方向弹性模量
            double Nu12 =compositeMaterialProperties[2]; // 单向带面内泊松比
            double G12 = compositeMaterialProperties[3];
            double Rho = compositeMaterialProperties[4];
            double alp1 = compositeMaterialProperties[5]; // 纵向膨胀系数
            double alp2 = compositeMaterialProperties[6]; // 横向膨胀系数

            double Nu21 = Nu12 * E2 / E1; // 单向带面内泊松比
            // 定义rad存放弧度，定义Q矩阵， 定义Z数组存放铺层Z坐标
            double[] rad = new double[layupIntArray.length];
            double[] Q11 = new double[layupIntArray.length];
            double[] Q12 = new double[layupIntArray.length];
            double[] Q22 = new double[layupIntArray.length];
            double[] Q16 = new double[layupIntArray.length];
            double[] Q26 = new double[layupIntArray.length];
            double[] Q66 = new double[layupIntArray.length];
            double[] Z = new double[layupIntArray.length+2];
            // 折减刚度矩阵
            double E11 = E1 / (1 - Nu12 * Nu21);
            double E22 = E2 / (1 - Nu12 * Nu21);
            double E12 = Nu21 * E1 / (1 - Nu12 * Nu21);
            double E66 = G12;
            double A11 = 0, A12 = 0, A22 = 0, A16 = 0, A26 = 0, A66 = 0;
            double B11 = 0, B12 = 0, B22 = 0, B16 = 0, B26 = 0, B66 = 0;
            double D11 = 0, D12 = 0, D22 = 0, D16 = 0, D26 = 0, D66 = 0;
            double thermoX = 0, thermoY = 0, thermoXY = 0;

            double thickness = Arrays.stream(layerThicknessArray).sum(); // 计算铺层总厚度
            Z[0] = -1 * thickness / 2; // Z0坐标值
            // 中间厚度所在层
            double currentThickness = 0;
            int currentLayer = 0;
            while (currentThickness < -Z[0]) {
                currentLayer++;
                currentThickness += layerThicknessArray[currentLayer];
            }
            // 处理不对称情况
            double upperThickness = currentThickness + Z[0];
            double lowerThickness = layerThicknessArray[currentLayer] - upperThickness;

            for (int i = layerThicknessArray.length-1; i >= currentLayer +2; i--){
                layerThicknessArray[i] = layerThicknessArray[i - 1];
            }
            layerThicknessArray[currentLayer] = lowerThickness;
            layerThicknessArray[currentLayer+1] = upperThickness;

            // 计算每一铺层的Z坐标
            for (int i = 1; i < layerThicknessArray.length; i++) {
                Z[i] = layerThicknessArray[i] + Z[i-1];
            }
            for (int i = currentLayer; i < layerThicknessArray.length; i++) {
                Z[i] = Z[i+1];
            }

            // 铺层角度、三角函数
            for (int i = 1; i < rad.length; i++) {
                rad[i] = layupIntArray[i] * Math.PI / 180;
                double sinRad = Math.sin(rad[i]);
                double cosRad = Math.cos(rad[i]);

                Q11[i] = E11 * Math.pow(cosRad, 4) +
                    2 * Math.pow(sinRad, 2) * Math.pow(cosRad, 2) * (E12 + 2 * E66) +
                    E22 * Math.pow(sinRad, 4);
                Q12[i] = (E11 + E22 - 4 * E66) * Math.pow(sinRad, 2) * Math.pow(cosRad, 2) +
                    E12 * (Math.pow(cosRad, 4) + Math.pow(sinRad, 4));
                Q22[i] = E11 * Math.pow(sinRad, 4) +
                    2 * Math.pow(sinRad, 2) * Math.pow(cosRad, 2) * (E12 + 2 * E66) +
                    E22 * Math.pow(cosRad, 4);
                Q16[i] = (E11 - E12 - 2 * E66) * sinRad * Math.pow(cosRad, 3) +
                    (E12 - E22 - 2 * E66) * cosRad * Math.pow(sinRad, 3);
                Q26[i] = (E11 - E12 - 2 * E66) * cosRad * Math.pow(sinRad, 3) +
                    (E12 - E22 - 2 * E66) * sinRad * Math.pow(cosRad, 3);
                Q66[i] = (E11 + E22 - 2 * E12 - 2 * E66) * Math.pow(sinRad, 2) * Math.pow(cosRad, 2) +
                    E66 * (Math.pow(cosRad, 4) + Math.pow(sinRad, 4));

                double Z1 = Z[i - 1];
                double Z2 = Z[i];

                A11 = A11 + Q11[i] * (Z2 - Z1);
                A12 = A12 + Q12[i] * (Z2 - Z1);
                A22 = A22 + Q22[i] * (Z2 - Z1);
                A16 = A16 + Q16[i] * (Z2 - Z1);
                A26 = A26 + Q26[i] * (Z2 - Z1);
                A66 = A66 + Q66[i] * (Z2 - Z1);

                B11 = B11 + Q11[i] * (Math.pow(Z2, 2) - Math.pow(Z1, 2)) / 2;
                B12 = B12 + Q12[i] * (Math.pow(Z2, 2) - Math.pow(Z1, 2)) / 2;
                B22 = B22 + Q22[i] * (Math.pow(Z2, 2) - Math.pow(Z1, 2)) / 2;
                B16 = B16 + Q16[i] * (Math.pow(Z2, 2) - Math.pow(Z1, 2)) / 2;
                B26 = B26 + Q26[i] * (Math.pow(Z2, 2) - Math.pow(Z1, 2)) / 2;
                B66 = B66 + Q66[i] * (Math.pow(Z2, 2) - Math.pow(Z1, 2)) / 2;

                D11 = D11 + Q11[i] * (Math.pow(Z2, 3) - Math.pow(Z1, 3)) / 3;
                D12 = D12 + Q12[i] * (Math.pow(Z2, 3) - Math.pow(Z1, 3)) / 3;
                D22 = D22 + Q22[i] * (Math.pow(Z2, 3) - Math.pow(Z1, 3)) / 3;
                D16 = D16 + Q16[i] * (Math.pow(Z2, 3) - Math.pow(Z1, 3)) / 3;
                D26 = D26 + Q26[i] * (Math.pow(Z2, 3) - Math.pow(Z1, 3)) / 3;
                D66 = D66 + Q66[i] * (Math.pow(Z2, 3) - Math.pow(Z1, 3)) / 3;

                thermoX = thermoX + (Z2 - Z1) * (Math.pow(cosRad, 2) * E11 * (alp1 + Nu21 * alp2) + Math.pow(sinRad, 2) * E22 * (Nu12 * alp1 + alp2));
                thermoY = thermoY + (Z2 - Z1) * (Math.pow(sinRad, 2) * E11 * (alp1 + Nu21 * alp2) + Math.pow(cosRad, 2) * E22 * (Nu12 * alp1 + alp2));
                thermoXY = thermoXY + (Z2 - Z1) * cosRad * sinRad * (E22 * (Nu12 * alp1 + alp2) - E11 * (alp1 + Nu21 * alp2));
                //thermoXY = 3.46944695195361E-17;
            }

            double[][] AData = {
                {A11, A12, A16},
                {A12, A22, A26},
                {A16, A26, A66}
            };
            RealMatrix A = MatrixUtils.createRealMatrix(AData);

            // 计算矩阵A的逆矩阵
            RealMatrix inverseMatrixA = new LUDecomposition(A).getSolver().getInverse();

            // 提取逆矩阵的元素
            double MA11 = inverseMatrixA.getEntry(0, 0);
            double MA22 = inverseMatrixA.getEntry(1, 1);
            double MA66 = inverseMatrixA.getEntry(2, 2);
            double MA12 = inverseMatrixA.getEntry(0, 1);

            // 计算整个铺层等效弹性模量 Exx_eqv, Eyy_eqv, Gxy_eqv 和泊松比 NU12_eqv, NU21_eqv
            double ExxEquivalent = 1 / MA11 / thickness;
            double EyyEquivalent = 1 / MA22 / thickness;
            double GxyEquivalent = 1 / MA66 / thickness;
            double NU12Equivalent = -1 * MA12 / MA11;
            double NU21Equivalent = -1 * MA12 / MA22;
            // 计算等效热膨胀系数
            double[][] thermoCoeData = {
                {thermoX, thermoY, thermoXY}
            };
            RealMatrix thermoCoe = MatrixUtils.createRealMatrix(thermoCoeData);

            RealMatrix thermoCoeTransposed = thermoCoe.transpose();
            RealMatrix matrixAlpEquivalent = inverseMatrixA.multiply(thermoCoeTransposed);

            double[] alpEquivalent = new double[3];
            alpEquivalent[0] = matrixAlpEquivalent.getEntry(0, 0);
            alpEquivalent[1] = matrixAlpEquivalent.getEntry(1, 0);
            alpEquivalent[2] = matrixAlpEquivalent.getEntry(2, 0);
            // 保存信息
            equivalentProperties[0] = ExxEquivalent;
            equivalentProperties[1] = EyyEquivalent;
            equivalentProperties[2] = GxyEquivalent;
            equivalentProperties[3] = thickness;
            equivalentProperties[4] = Rho;
            equivalentProperties[5] = alpEquivalent[0];
            equivalentProperties[6] = alpEquivalent[1];
            equivalentProperties[7] = alpEquivalent[2];

        }
    }

    private static String[] expandSymmetrically(String[] originalArray, int symmetrical){
        if (symmetrical == 0) {
            return originalArray;
        }
        int length = originalArray.length;
        String[] symmetricArray = new String[length*2];
        System.arraycopy(originalArray, 0, symmetricArray, 0, length);
        for (int i = 0; i < length; i++) {
            symmetricArray[length+i] = originalArray[length-1-i];
        }
        return symmetricArray;
    }

    public static void main(String[] args) {

        String plyCode = "1423";
        double singleThickness = 0.187;

        int symmetricalA = 0;
        String layupSequenceA = "45,0,0,-45,90,-45,0,0,45";

        int symmetricalB = 1;
        String layupSequenceB = "45,0,0,-45,0,0,90";

        int symmetricalC = 0;
        String layupSequenceC = "-45,45,45,90,-45,0,0,-45,90,45,45,-45";

        int symmetricalD = 0;
        String layupSequenceD = "45,-45,0,45,0,0,-45,0,0,90,90,0,0,-45,0,0,45,0,-45,45";

        int symmetricalE = 0;
        String layupSequenceE = "45,-45,0,45,0,-45,0,90,90,0,-45,0,45,0,-45,45";

        double[] compositeMaterialProperties = {157000, 9650, 0.31, 4600,
                                                0.00000000154, -0.0000016, 0.0000513};
        double[] equivalentProperties = new double[8];

        calculateEquivalentMaterial(plyCode, singleThickness, symmetricalA, layupSequenceA,
            compositeMaterialProperties, equivalentProperties);
        System.out.println("序列号923: " + Arrays.toString(equivalentProperties));

        calculateEquivalentMaterial(plyCode, singleThickness, symmetricalB, layupSequenceB,
            compositeMaterialProperties, equivalentProperties);
        System.out.println("序列号1423: " + Arrays.toString(equivalentProperties));

        calculateEquivalentMaterial(plyCode, singleThickness, symmetricalC, layupSequenceC,
            compositeMaterialProperties, equivalentProperties);
        System.out.println("序列号1221: " + Arrays.toString(equivalentProperties));

        calculateEquivalentMaterial(plyCode, singleThickness, symmetricalD, layupSequenceD,
            compositeMaterialProperties, equivalentProperties);
        System.out.println("序列号2023: " + Arrays.toString(equivalentProperties));

        calculateEquivalentMaterial(plyCode, singleThickness, symmetricalE, layupSequenceE,
            compositeMaterialProperties, equivalentProperties);
        System.out.println("序列号1623: " + Arrays.toString(equivalentProperties));


    }
}
